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Abstract. Weak Galerkin (WG) refers to general finite element methods for partial differential 
equations in which differential operators are approximated by weak forms through the usual inte- 
gration by parts. In particular, WG methods allow the use of discontinuous finite element functions 
in the algorithm design. One of such examples was recently introduced in [54] for solving second 
order elliptic problems. The goal of this paper is to apply the WG method of [54] to the Helmholtz 
equation with high wave numbers. Several test scenarios are designed for a numerical investigation 
on the accuracy, convergence, and robustness of the WG method in both inhomogeneous and ho- 
mogeneous media over convex and non-convex domains. Our numerical experiments indicate that 
weak Galerkin is a finite element technique that is easy to implement, and provides very accurate 
and robust numerical solutions for the Helmholtz problem with high wave numbers. 

Key words. Galerkin finite element methods, discrete gradient, the Helmholtz equation, weak 
Galerkin 
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1. Introduction. In this paper, we explore the use of a weak Galerkin (WG) 
finite element method for solving the nonhomogeneous Helmholtz equation with high 
wave numbers 



where k is the wave number, / represents a harmonic source, and d = d{x^ y) is a spa- 
tial function describing the dielectric properties of the medium. The Helmholtz equa- 
tion governs many macroscopic wave phenomena in the frequency domain in- 
cluding wave propagation, guiding, radiation and scattering, where the time-harmonic 
behavior can be assumed. The numerical solution to the Helmholtz equation plays a 
vital role in a wide range of applications in electromagnetics, optics, and acoustics, 
such as antenna analysis and synthesis, radar cross section calculation, simulation of 
ground or surface penetrating radar, design of optoelectronic devices, acoustic noise 
control, and seismic wave propagation. However, it remains a challenge to design 
robust and efficient numerical algorithms for the Helmholtz equation, especially when 
high wave numbers or highly oscillatory solutions are involved [58 . 

Physically, the Helmholtz problem is usually defined on an unbounded exterior 
domain with the so-called Sommerfeld radiation condition holding at infinity [381 153] 
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where i = is the imaginary unit and r is the radial direction. Here we have 

assumed that d = 1 in far field. In computational electromagnetics, the exterior 
domain problems are often solved numerically by introducing a bounded domain Q 
with an artificial boundary and imposing certain boundary conditions on so 
that nonphysical reflections from the boundary can be eliminated or minimized. For 
the Helmholtz exterior problems, the non-reflecting condition is commonly chosen 
as a Dirichlet-to-Neumann (DtN) mapping, which relates the wave solution to its 
derivatives on dQ |38l |53| 



(1.3) d\/u • n - T{u) = g, on dQ, 

where n denotes the outward normal direction of dfl^ T is a DtN integral operator, and 
^ is a given data function. For sufliciently large the nonlocal boundary condition 
( |1.3[ ) can be approximated by a Robin boundary condition [38[ |53] 

(1.4) d\/u • n — iku = on dQ, 

which is essentially a flrst order absorbing boundary condition. 

In the present paper, we consider the following prototype Helmholtz problem 

(1.5) - V • (dVu) -k^u = f in 

(1.6) d\/u • n — iku = g on dQ. 

Finite element methods for such Helmholtz problems can be classifled as two cate- 
gories. The flrst category consists of methods that use continuous functions to approx- 
imate the solution u and the other refers to methods with discontinuous approximation 
functions. 

Continuous Galerkin (CG) flnite element methods employ continuous piecewise 
polynomials to approximate the true solution of (1.5)-(1.6) and lead to a simple for- 
mulation: flnd Uh ^ Vh C H^{Q) satisfying 

(1.7) {dVuh, Vvh) - k^{uh, Vh) + ik{uh, Vh)dn = (/, Vh) + (^, Vh)dn 

for all Vh G Vh^ where Vh is an properly deflned flnite element space consisting of 
continuous piecewise polynomials. 

Discontinuous Galerkin (DG) flnite element methods for the Helmholtz equation 
with d = 1 (e.g. [33j) seek Uh ^ Vh C L'^{Q) satisfying 

Y^{VUh, VVh)T- Y^{{{^Uh}, [Vh])e + ^({VV^}, [Uh])e) - k\uh, Vh) 
T e 

Mk{Uh, Vh)dn^iC^^lK^{[^h], b/i])e + ^<^2^e([^^], I^De 



(1.8) +E^3/^e-'([^], [^])e) = U.Vh) + {g.Vh)8n. 



for all Vh ^Vh^ where a^, z=l,2,3 are penalty parameters and Vh is a space of discon- 
tinuous piecewise polynomials. 

The continuous Galerkin flnite element formulation (1.7) is natural and easy to 
implement. However, the errors of continuous Galerkin flnite element solutions de- 
teriorate rapidly when k becomes large. On the other hand, the formulation of DG 
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methods (1.8) is complex, and the stabihty and accuracy of DG methods heavily de- 
pend on the selection of penalty parameters (a good determination for the parameter 
values has been an issue for DG methods). 

The objective of the present paper is to introduce a weak Galerkin (WG) finite 
element method for solving the Helmholtz problem with high wave numbers and to 
investigate the robustness and effectiveness of such WG method through many care- 
fully designed numerical experiments. The weak Galerkin finite element formulation 
was first developed in [54] for solving the second order elliptic equations. Through 
rigorous error analysis, optimal order of convergence of the WG solution in both dis- 
crete norm and norm is established under minimum regularity assumptions in 
[54] , Nevertheless, no experimental results have been reported in [54] . 

Generally speaking, like the DG method, the WG finite element method allows one 
to use discontinuous functions in the finite element procedure. The main idea behind 
weak Galerkin methods lies in the approximation of differential operators by weak 
forms for discontinuous finite element functions defined on a partition of the domain. 
For example, for the gradient operator, one may introduce a weak gradient operator 
\/d over discontinuous functions Vh = {vo^Vb}, where '^o is defined in the interior of 
the element and v^j is defined on the boundary of the element. The weak gradient 
operator will then be employed to form a weak Galerkin finite element formulation 



for ( |1.5| )-( |L6l ): find Uh G Vh such that for all Vh G Vh we have 

(1.9) {dVdUh, ^dVh) - k^{uo, vo) + ik{ub, Vb)QQ = (/, vq) + {g, Vb)Qn. 



The discrete gradient \/dV in ( |l.9[ ) is a variable locally calculated on each element. 
The weak Galerkin methods have many nice features. First, the formulation of WG 
method is simple, easy to implement, and involves no penalty parameters for users 
to select. The weak Galerkin is equipped with elements of polynomials of any degree 
k > 0. Secondly, the weak Galerkin method conserves mass with a well defined 
numerical fiux. In some sense, weak Galerkin finite element methods enjoy both the 
simplicity of CG method and the flexibility of DG method. 

In the present numerical investigation, we are particularly interested in the perfor- 
mance of the WG method for solving the Helmholtz equation with high wave numbers. 
In general, the numerical performance of any finite element solution to the Helmholtz 
equation depends significantly on the wave number k. When k is very large - repre- 
senting a highly oscillatory wave, the mesh size h has to be sufliciently small for the 
scheme to resolve the oscillations. To keep a fixed grid resolution, a natural rule is to 
choose kh to be a constant in the mesh refinement, as the wave number k increases 
[4T1IT2]. However, it is known [4TJ|42j|8ll3 that, even under such a mesh refinement, 
the errors of continuous Galerkin finite element solutions deteriorate rapidly when k 
becomes larger. This non-robust behavior with respect to k is known as the "pollution 
effect" [4l]|42] ISl E] . Usually, under the small magnitude assumption of kh, the rela- 
tive error of a pih order finite element solution in the i^^ -semi- norm consists of two 
parts [4TJ|42], i.e., an error term of the best approximation behavior like 0(/c^/i^) and 
a pollution error term behavior like 0{pP^^h'^P). The error of best approximation 
is essentially due to the interpolation error on a discretized grid and is of bounded 
magnitude if kh = constant. Nevertheless, the pollution error term dominates when 
k is large and is responsible for the non-robustness behavior in the finite element 
solutions to the Helmholtz equation. 

To the end of eliminating or substantially reducing the pollution errors, various 
numerical approaches have been developed in the literature for solving the Helmholtz 
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equation. For the continuous and least-squares finite element methods, a popular 
strategy [HI |9l [45l [TOl |46] to reduce the pollution error is to include some analytical 
knowledge of the problem, such as characteristics of asymptotic or exact solution, in 
the finite element space. The resulting local basis functions of non-polynomial shape 
yield improved performance [9l |45l [TOl [46]. Similarly, analytical information is 
incorporated in the basis functions of the boundary element methods to address the 
high frequency problems [36 , 44 , 23 . Based on the geometrical optics and geometrical 
theory of diffraction, asymptotic solutions of the Helmholtz equation that represent 
important wave propagation directions are coupled with high frequency oscillations to 
built boundary element approximation spaces. Consequently, the number of degrees 
of freedom of such boundary element methods virtually does not depend on the wave 
number k for the Helmholtz equation [36l [44l [23]. In [22 , plane wave solutions 
traveling in a large number of directions have also been employed as basis functions in 
the ultra weak variational formulation (UWVF) [21] . This allows the use of a coarse 
mesh when resolving high oscillatory wave solutions. The UWVF method can be 
regarded as an upwinding DG discretization for a first order system obtained through 
the introduction of an adjoint field. Based on such a viewpoint, error estimates of 
the UWVF on the entire domain have been recently presented in [T9^ . Multiple plane 
wave basis functions are also employed in a DG method using both low order elements 
[3T] and high order elements [32] . The weak continuity of the solution at the element 
boundaries is enforced by introducing a Lagrange multiplier. Being more robust than 
classical variational approaches, this DG method can resolve high frequency short 
wave problems very well with a ratio of up to three wavelengths per element [31] |32] . 

Besides changing basis functions, there are other improvements that can be made 
in the DG framework to reduce the pollution error. By using piecewise linear polyno- 
mials as the basis functions, stable interior penalty DG methods have been developed 
in [33] through penalizing not only the jumps of the function values, but also the 
jumps of the normal and tangential derivatives across the element edges. Robust 
results against pollution effect can be achieved through a careful selection of penalty 
parameters [33]. Also targeting on enhancing stability, DG methods can be estab- 
lished by formulating the wave equation as a first-order system [6[ |25] . The classical 
techniques, such as the enforcement of weak continuity via fiuxes form, can be used in 
such DG solutions to the Helmholtz equation. The dispersive and dissipative behavior 
of selected DG schemes for both second order Helmholtz equation and the correspond- 
ing first order system have been examined in [2] . Hybridizable DG method is another 
efficient finite element method for the Helmholtz equation [29] |30] . By introducing 
new degree of freedoms on the boundary of elements, parametrization is conducted 
element by element so that the final linear system consists of unknowns only from the 
skeleton of the mesh. This greatly reduces the size of linear systems compared to the 
standard DG scheme. The error analysis [37] shows that the condition number of the 
condensed matrix of the hybridizable DG methods [29l [30] could be independent of 
the wave number. 

The pollution effect can also be alleviated by controlling the numerical dispersion, 
i.e., the phase difference between the numerical and exact waves. This is because the 
pollution error is known to be directly related to the dispersion error [40l [1] , while 
the best approximation error is non-dispersive. Thus, the reduction of the dispersive 
error is equivalent to the reduction of the pollution error. Consequently, high order 
methods, such as local spectral method [12 , spectral Galerkin methods [52l[53], and 
spectral element methods [39l |3] are less vulnerable to the pollution effect, because 
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they produce negligible dispersive errors. 

The rest of this paper is organized as follows. In Section 2, we will introduce 
a weak Galerkin finite element method for the Helmholtz equation by following the 
idea presented in [54 . In Section 3, we will present some error estimates for the 
WG method for the Helmholtz equation. Finally in Section 4, we shall present some 
numerical results obtained from the weak Galerkin method with various orders. 

2. A Weak Galerkin Finite Element Method. Let 7/i be a partition of the 
domain Q with mesh size h. Assume that the partition Th is shape regular so that the 
routine inverse inequality in the finite element analysis holds true (see [26 ). Define 
{v,w)k = Jj^vwdx and {v^w)qk = J^j^vwds. 

For each triangle T G 7/i, let and dT denote the interior and boundary of 
T respectively. Denote by Pj{T^) the set of polynomials in with degree no more 
than j, and Pi{dT) the set of polynomials on each segment (edge or face) of dT with 
degree no more than I. We emphasize that functions of Pi{dT) are defined on each 
edge/face and there is no continuity required across different edges/faces. Define a 
global function v = {vo^Vb} where '^o and represent the values of v on T^ and 
dT respectively for each T e Th- Now we are ready to define a weak Galerkin finite 
element space as follows 

(2.1) Vh = {v = {vo,Vb} : {vo,Vb}\T G Pj{T'') x P,(5T),VT G %} . 

For each v = {vq^ v^} G Vh^ we define the discrete gradient of v on each element T by 
the following equation: 



(2.2) 



/ VdV-qdT = - / voV-qdT^ [ Vbq-nds, yqeVr{T), 

JT JT JdT 



where VriT) is a subspace of the set of vector- valued polynomials of degree no more 
than r on T. 

For the purpose of easy demonstration, we consider a special 2-D weak Galerkin 
element of Vn and Vr{T) with j = l = k>{)m and Vr{T) = RTk{T). More weak 
Galerkin elements can be found in [54]. Here RTk(T) is the usual Raviart-Thomas 
element [47 of order k which has the form 

Rn{T) = PkiTf ^ PkiT)^, 

where Pk{T) is the set of homogeneous polynomials of degree k and x = (xi, X2). 

Weak Galerkin Algorithm 1. A numerical approximation for l[l.^ -{1.6) 
can be obtained by seeking = {uq^ u^} G Vh such that for all = {'^o, ^b} ^ Vh 

(2.3) {dVdUh, ^d^h) - k^{uQ, vq) + ik{ub, Vb)dn = (/, ^o) + (^, Vb)dn' 



In the following, we will use the lowest order weak Galerkin element {k=0) as 
an example to demonstrate how one might implement weak Galerkin finite element 
method for solving the Helmholtz problem (1.5)-(1.6). Let N{T) and N{e) denote the 
number of triangles and the number of edges associated with triangulation Th- Let 
Sh denote the union of the boundaries of the triangles T of Tk- Let be a function 
which takes value one in the interior of triangle and zero everywhere else. Let tpj 
be a function that takes value one on edge ej and zero everywhere else. The weak 
Galerkin finite element space Vh with k = has the form 



Vh = span{(/)i,--- ,(/)Ar(T),V^i,--- ,V^Ar(e)}- 
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The methodology of implementing WG is the same as that for continuous Galerkin 
finite element method except that the standard gradient operator V must be replaced 
by the discrete gradient operator V^^. The discrete gradient operator is an exten- 
sion of the standard gradient operator for smooth functions to non-smooth functions. 
Therefore, the key step of implementing weak Galerkin method is to compute a local 
variable V^'^ for v = {vo^Vb} G Vh element by element. For the case /c = 0, we have 
\/dV G RTo{T) on each element T G Th where 

RTo{T)=( ItZ) = spaniel, ^1,^3}. 



b -\- cy 

For example, we can choose Oi as follows 



Thus on each element T G 7/i, V^^'^ = Yl^^iCjOj. Using the definition of discrete 



gradient (2.2), we find Cj by solving the following linear system: 



{0i,Oi)t {e2,0i)T {03,0i)t \ / ci \ / -{vo,V ■ 0i)t + {vb, ^1 • n)aT 

i01,O2)T {02,02)t {03,02)t ^2 = -(uq, V • ^2)t + (^6, ^2 ' n)aT 

(^1,^3)t {02,03)t {03,93)t J \ C3 J V "(vo, V • ^3)t + (^6, ^3 ' n)aT 

The inverse of the above coefficient matrix can be obtained explicitly or numerically 
through a local matrix solver. Once S/d^^i S/d'ipj are computed, the rest of steps 
are the same as those for the classical finite element method. 

3. Error Estimates. Denote by QhU = {Qqu^ Qb^} the projection onto 
Pj{T^) X Pj{dT). In other words, on each element T, the function Qqu is defined as 
the projection of u in Pj{T) and on 9T, Qi^u is the projection in Pj{dT). 

Consider the following Helmholtz problem: 

(3.1) - V • (dVu) -k^u = f in ft, 

(3.2) u = 9 on dQ. 



For weak Galerkin finite element space Vh defined in (2.1 ), we define a subspace 
of Vh with vanishing boundary values on dfl; i.e., 

Vh '= {v = {vo.Vb} e Vh,vi,\9TndQ = 0, VT G %} • 



The weak Galerkin finite element method for the Helmholtz problem (3.1)-(3.2) is 
stated as follows. 

Weak Galerkin Algorithm 2. A numerical approximation for ^3J\ ) and p.^P 
can be obtained by seeking Uh = {uq^ui)} G Vh satisfying U}) = Q^g on and the 
following equation: 

(3.3) {dVdUh, Vdv) - k^{uQ, vq) = (/, vq), W v = {vo.Vb} G Vh, 



For a sufficiently small mesh size h, the following error estimate holds true. A 
proof of this error estimate can be found in [54] . 



Theorem 3.1. Let u G H^{Vt) be the solution (3.1) and (3.2), and Uh be a weak 
Galerkin approximation of u arising from {3.^. Assume that the exact solution u is 
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sufficiently smooth such that u G H'^^^{Q) with < m < /c + 1. Then, there exists a 
constant C such that 

(3.4) \\Vd{uh - Qhu)\\ < C(/i™||u|U+i + h^+-'\\f - Qo/ll) 

(3.5) \\uh - Qhu\\ < C(/i™+^||m|U+i + h'+'\\f - Qof\\), 

provided that the dual problem of (3.1)- (3.2) has the H^^^(Q) regularity with s G 



(0,1]. 

4. Numerical Experiments. In this section, we examine the WG method 
by testing its accuracy, convergence, and robustness for solving two dimensional 
Helmholtz equations. The pollution effect due to large wave numbers will be par- 
ticularly investigated and tested numerically. For convergence tests, both piecewise 
constant and piecewise linear finite elements will be considered. To demonstrate WG's 
robustness, the Helmholtz equation in both homogeneous and inhomogeneous media 
will be solved on convex and non-convex computational domains. For simplicity, a 
structured triangular mesh is employed in all cases. The mesh generation and all 
computations are conducted in the MATLAB environment. 

Two types of relative errors are measured in our numerical experiments. The first 
one is the relative LF' error and the second is the relative error. The relative 
error is defined by 

\\uh - Qhu\\ 

WQhuW ' 

where is the WG approximation to the solution and QhU = Qh"^} is the 

projection onto Pj{T^) x Pj{dT). In other words, on each element T, the function 
Qqu is defined as the projection of u in PjiT) and on 9T, Qi)U is the LP' projection 
in Pj{dT). The relative error is defined in terms of the discrete gradient 

\\Vd{uh-Qhu)\\ 

WdQhuW 

However, in the present study, the i^^-semi-norm will be calculated as 

\\\uh - Qhu\f = h~^{uo -Ub- {QoU - QbU),Uo - Uj, - {QqU - QhU))e 

for the lowest order finite element (i.e., piecewise constants). For piecewise linear 
elements, we use the original definition of V^^ to compute the -^-semi-norm \\V d{uh — 

Qhu)l 

4.1. A convex Helmholtz problem. We first consider a homogeneous Helmholtz 
equation defined on a convex hexagon domain, which has been studied in [33 . The 
domain Vt is the unit regular hexagon domain centered at the origin (0,0), see Fig. 



4.1 



(left). Here we set d = 1 and / = sin(/cr)/r in (1.5), where r = y^x^+^- The 
boundary data g in the Robin boundary condition (1.6) is chosen so that the exact 
solution is given by 

. cos( kr) cos k-\-i sink ^ . 

where J^{z) are Bessel functions of the first kind. Let Th denote the regular triangu 



lation that consists of triangles of size = l/A/", as shown in Fig. 4.1 (left) for 
Ti. 



Fig. 4.1. Geometry of testing domains and sample meshes. Left: a convex hexagon domain; 
Right: a non- convex imperfect circular domain. 



Table 4.1 

Convergence of piecewise constant WG for the Helmholtz equation on a convex domain with 
wave number k = 1. 



h 


relative 


relative 


error 


order 


error 


order 


5.00e-01 


2.49e-02 




4.17e-03 




2.50e-01 


l.lle-02 


1.16 


1.05e-03 


1.99 


1.25e-01 


5.38e-03 


1.05 


2.63e-04 


2.00 


6.25e-02 


2.67e-03 


1.01 


6.58e-05 


2.00 


3.13e-02 


1.33e-03 


1.00 


1.64e-05 


2.00 


1.56e-02 


6.65e-04 


1.00 


4.11e-06 


2.00 



Table 4.2 

Convergence of piecewise linear WG for the Helmholtz equation on a convex domain with wave 
number k = 5. 



h 


relative 


relative 


error 


order 


error 


order 


2.50e-01 


9.48e-03 




2.58e-04 




1.25e-01 


2.31e-03 


2.04 


3.46e-05 


2.90 


6.25e-02 


5.74e-04 


2.01 


4.47e-06 


2.95 


3.13e-02 


1.43e-04 


2.00 


5.64e-07 


2.99 


1.56e-02 


3.58e-05 


2.00 


7.06e-08 


3.00 


7.81e-03 


8.96e-06 


2.00 


8.79e-09 


3.01 



Table 4.1 illustrates the performance of the WG method with piecewise constant 
elements for the Helmholtz equation with wave number k = 1. Uniform triangular 
partitions were used in the computation through successive mesh refinements. The 
relative errors in norm and semi-norm can be seen in Table l4.ll The Table 
also includes numerical estimates for the rate of convergence in each metric. It can 
be seen that the order of convergence in the relative semi-norm and relative 
norm are, respectively, one and two for piecewise constant elements. 



0.5 



-0.5 



-1 




■0.1 -0-5 



-0.3 _i 




0.5 0.5 1 



Fig. 4.2. WG solutions for the non-convex Helmholtz problem with k = 4 and ^ = 1. Left: 
Mesh level 1; Right: Mesh level 6. 



High order of convergence can be achieved by using corresponding high order 
finite elements in the present WG framework. To demonstrate this phenomena, we 
consider the same Helmholtz problem with a slightly larger wave number k = 5. The 
WG with piecewise linear functions was employed in the numerical approximation. 
The computational results are reported in Table |4.2[ It is clear that the numerical 
experiment validates the theoretical estimates. More precisely, the rates of conver- 
gence in the relative semi-norm and relative norm are given by two and three, 
respectively. 

4.2. A non-convex Helmholtz problem. We next explore the use of the WG 

method for solving a Helmholtz problem defined on a non-convex domain, see Fig. 
|4.1| (right). The medium is still assumed to be homogeneous, i.e., d = 1 in (1.5). We 



are particularly interested in the performance of the WG method for dealing with the 
possible field singularity at the origin. For simplicity, only the piecewise constant RTq 
elements are tested for the present problem. Following [37 , we take / = in (1.5) 
and the boundary condition is simply taken as a Dirichlet one: = ^ on dfl. Here g 
is prescribed according to the exact solution [3^ 



(4.2) 



y'^) cos(^ arctan(?//x)). 



In the present study, the wave number was chosen as /c = 4 and three values 
for the parameter ^ are considered; i.e., ^ = 1,^ = 3/2 and ^ = 2/3. The same 
triangular mesh is used in the WG method for all three cases. In particular, an initial 
mesh is first generated by using MATLAB with default settings, see Fig. |4.1| (right). 
Next, the mesh is refined uniformly for five times. The WG solutions on mesh level 



1 and mesh level 6 are shown in Fig. |4.2[ Fig. |4.3[ and Fig. |4.4[ respectively, for 
^ = 1, ^ = 3/2, and ^ = 2/3. Since the numerical errors are quite small for the 
WG approximation corresponding to mesh level 6, the field modes generated by the 
densest mesh are visually indistinguishable from the analytical ones. In other words. 



the results shown in the right charts of Fig. |4.2[ Fig. |4.3[ and Fig. |4.4| can be regarded 
as analytical results. It can be seen that in all three cases, the WG solutions already 
agree with the analytical ones at the coarsest level. Moreover, based on the coarsest 
mesh, the constant function values can be clearly seen in each triangle, due to the use 
of piecewise constant RTq elements. Nevertheless, after the initial mesh is refined for 
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Fig. 4.3. WG solutions for the non-convex Helmholtz problem with k = 4 and ^ = 3/2. Left: 
Mesh level 1; Right: Mesh level 6. 
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Fig. 4.4. WG solutions for the non-convex Helmholtz problem with k = 4 and ^ = 2/3. Left: 
Mesh level 1; Right: Mesh level 6. 



Table 4.3 

Numerical convergence test for the non-convex Helmholtz problem with k = 4 and ^ = 1. 



h 


relative 


relative 


error 


order 


error 


order 


2.44e-01 


5.64e-02 




1.37e-02 




1.22e-01 


2.83e-02 


1.00 


3.56e-03 


1.95 


6.10e-02 


1.42e-02 


0.99 


8.98e-04 


1.99 


3.05e-02 


7.14e-03 


1.00 


2.25e-04 


2.00 


1.53e-02 


3.57e-03 


1.00 


5.63e-05 


2.00 


7.63e-03 


1.79e-03 


1.00 


1.41e-05 


2.00 



five times, the numerical plots shown in the right charts are very smooth. A perfect 
symmetry with respect to the x-axis is clearly seen. 

We next investigate the numerical convergence rates for WG. The numerical errors 
of the WG solutions for ^ = 1, ^ = 3/2 and ^ = 2/3 are listed, respectively, in Table 
4.3, Table 4.4, and Table 4.5 It can be seen that for ^ = 1 and ^ = 3/2, the 
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Table 4.4 

Numerical convergence test for the non-convex Helmholtz problem with k = 4 and ^ = 3/2. 



h 


relative 


relative 


error 


order 


error 


order 


2.44e-01 


5.56e-02 




1.12e-2 




1.22e-01 


2.81e-02 


0.98 


3.02e-03 


1.89 


6.10e-02 


1.42e-02 


0.99 


8.06e-04 


1.91 


3.05e-02 


7.14e-03 


0.99 


2.12e-04 


1.92 


1.53e-02 


3.58e-03 


1.00 


5.54e-05 


1.94 


7.63e-03 


1.79e-03 


1.00 


1.44e-05 


1.95 



Table 4.5 

Numerical convergence test for the non-convex Helmholtz problem with k = 4 and ^ = 2/3. 



h 


relative 


relative 


error 


order 


error 


order 


2.44e-01 


1.07e-01 




5.24e-02 




1.22e-01 


5.74e-02 


0.90 


2.18e-02 


1.27 


6.10e-02 


3.23e-02 


0.83 


9.01e-03 


1.27 


3.05e-02 


1.89e-02 


0.77 


3.68e-03 


1.29 


1.53e-02 


1.14e-02 


0.73 


1.49e-03 


1.31 


7.63e-03 


6.99e-03 


0.71 


5.96e-04 


1.32 



numerical convergence rates in the relative and errors remain to be first and 
second order, while the convergence orders degrade for the non-smooth case ^ = 2/3. 



Mathematically, for both ^ = 3/2 and ^ = 2/3, the exact solutions (4.2) are known to 
be non-smooth across the negative x-axis if the domain was chosen to be the entire 
circle. However, the present domain excludes the negative x-axis. Thus, the source 



term / of the Helmholtz equation (1.5) can be simply defined as zero throughout Q. 
Nevertheless, there still exists some singularities at the origin (0,0). In particular, it 
is remarked in [37] that the singularity lies in the derivatives of the exact solution at 
(0, 0). Due to such singularities, the convergence rates of high order DG methods are 
also reduced for ^ = 3/2 and ^ = 2/3 [37 . In the present study, we further note that 
there exists a subtle difference between two cases ^ = 3/2 and ^ = 2/3 at the origin. 
To see this, we neglect the second cos(-) term in the exact solution (|4.2|) and plot the 



Bessel function of the first kind J^{k\r\) along the radial direction r, see Fig. 4.5 It is 
observed that the Bessel function of the first kind is non-smooth for the case ^ = 2/3, 
while it looks smooth across the origin for the case ^ = 3/2. Thus, it seems that the 
first derivative of Js/2{k\r\) is still continuous along the radial direction. This perhaps 
explains why the present WG method does not experience any order reduction for the 
case ^ = 3/2. In [37 , locally refined meshes were employed to resolve the singularity 
at the origin so that the convergence rate for the case ^ = 2/3 can be improved. 
We note that local refinements can also be adopted in the WG method for a better 
convergence rate. A study of WG with grid local refinement is left to interested parties 
for future research. 

4.3. A Helmholtz problem with inhomogeneous media. We consider a 
Helmholtz problem with inhomogeneous media defined on a circular domain with 



radius R. Note that the spatial function d{x^y) in the Helmholtz equation (1.5) 
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Fig. 4.5. The Bessel function of the first kind J^(k\r\) across the origin. 



represents the dielectric properties of the underlying media. In particular, we have 
d = - in the electromagnetic applications [56 , where e is the electric permittivity. In 
the present study, we construct a smooth varying dielectric profile: 



(4.3) 



d(r) = -S{r) + -{l-S{r)), 



where r = y^x^ + ?/^, ei and 62 are dielectric constants, and 



(4.4) 



S{r) = { 



1 







if r < a. 



(^) +3(^) iia<r<b, 
if r > 6, 



with a < b < R. An example plot of d{r) and S{r) is shown in Fig. 4.6 In classical 
electromagnetic simulations, e is usually taken as a piecewise constant, so that some 
sophisticated numerical treatments have to be conducted near the material interfaces 
to secure the overall accuracy [56 . Such a procedure can be bypassed if one considers 
a smeared dielectric profile, such as (4.3). We note that under the limit 6 ^ a, 
a piecewise constant profile is recovered in (4.3). In general, the smeared profile 



(4.3) might be generated via numerical filtering, such as the so-called e-smoothing 
technique [51 in computational electromagnetics. On the other hand, we note that 
the dielectric profile might be defined to be smooth in certain applications. For 
example, in studying the solute-solvent interactions of electrostatic analysis, some 
mathematical models ^ have been proposed to treat the boundary between the 
protein and its surrounding aqueous environment to be a smoothly varying one. In 



fact, the definition of (4.3) is inspired by a similar model in that field [24 . 



In the present study, we choose the source of the Helmholtz equation (1.5) to be 



(4.5) 



/(r) = P[d{r) - 1] Jo(fcr) + kd'(r)Ji(kr) 



where 
(4.6) 



d'{r) = ( ^ - ) S'{r) 
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Fig. 4.6. An example plot of smooth dielectric profile d(r) and S(r) with a = 1, h = ?f and 
R = 5. The dielectric coefficients of protein and water are used, i.e., ei = 2 and 62 = 80. 



and 



(4.7) 



S'{r) 



if r < a, 







if 6 < r, 



For simplicity, a Dirichlet boundary condition is imposed at r 
g is prescribed according to the exact solution 



R with u = g. Here 



(4.8) 



u — Jo{kr). 



Our numerical investigation assumes the value of a = 1, 6 = 3 and R = 5. 
The wave number is set to be /c = 2. The dielectric coefficients are chosen as ei = 
2 and €2 = 80, which represents the dielectric constant of protein and water [24l 
|57] , respectively. The WG method with piecewise constant finite element functions 
is employed to solve the present problem with inhomogeneous media in Cartesian 



coordinate. Table 4.6 illustrates the computational errors and some numerical rate of 



convergence. It can be seen that the numerical convergence in the relative L error 
is not uniform, while the relative error still converges uniformly in first order. 
This phenomena might be related to the non-uniformity and smallness of the media 
in part of the computational domain. Nevertheless, the averaged convergence rate in 
the relative norm is about 2.12. Overall, we are confident that the WG method is 
accurate and robust in solving the Helmholtz equations with inhomogeneous media. 

4.4. Large wave numbers. We finally investigate the performance of the WG 
method for the Helmholtz equation with large wave numbers. The homogeneous 



Helmholtz problem studied in the Subsection 4.1 is employed again for this purpose. 
Also, the RTq and RTi elements are used to solve the homogeneous Helmholtz equa- 
tion with the Robin boundary condition. Since this problem is defined on a structured 
hexagon domain, a uniform triangular mesh with a constant mesh size h throughout 
the domain is used. This enables us to precisely evaluate the impact of the mesh 
refinements. Following the literature works [I2l|33], we will focus only on the relative 
semi-norm in the present study. 
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Table 4.6 



Numerical convergence test of the Helmholtz equation with inhomogeneous media 



h 


relative 


relative 


error 


order 


error 


order 


1.51e-00 


2.20e-01 




1.04e-00 




7.54e-01 


1.24e-01 


0.83 


1.20e-01 


3.11 


3.77e-01 


6.24e-02 


0.99 


1.81e-02 


2.73 


1.88e-01 


3.13e-02 


1.00 


5.71e-03 


1.67 


9.42e-02 


1.56e-02 


1.00 


2.14e-03 


1.42 


4.71e-02 


7.82e-03 


1.00 


5.11e-04 


2.06 



10 



10 



> -2 
10 'r 



10"^ 



10" 
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Fig. 4.7. Relative error of the WG solution. Left: with respect to 1/h; Right: with respect 
to wave number k. 



To study the non-robustness behavior with respect to the wave number /c, i.e., 
the pollution effect, we solve the corresponding Helmholtz equation by using piecewise 
constant WG method with various mesh sizes for four wave numbers /c = 5, /c = 10, 
k = 50, and k = 100, see Fig. |4.7| (left) for the WG performance. From Fig. 4.7 



(left), it can be seen that when h is smaller, the WG method immediately begins to 
converge for the cases k = b and k = 10. However, for large wave numbers /c = 50 
and k = 100, the relative error remains to be about 100%, until h becomes to be quite 
small or 1/h is large. This indicates the presence of the pollution effect. In the same 
figure, we also show the errors of different k values by fixing kh = 0.25. Surprisingly, 
we found that the relative error does not evidently increase as k becomes larger. 
The convergence line for kh = 0.25 looks almost flat, with a very little slope. We 
note that the present result of the WG method is as good as the one reported in 
^33j by using a penalized discontinuous Galerkin approach with optimized parameter 
values [33]. We emphasize that WG has advantage over the penalized DG in that no 
parameters are involved in the numerical scheme. 

On the other hand, the good performance of the WG method for the case kh = 
0.25 does not mean that the WG method could be free of pollution effect. In fact, 
it is known theoretically [9 that the pollution error cannot be eliminated completely 
in two- and higher-dimensional spaces for Galerkin finite element methods. In the 
right chart of Fig. |4.7[ we examine the numerical errors by increasing /c, under the 
constraint that kh is a constant. Huge wave numbers, up to /c = 240, are tested. 
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It can be seen that when the constant changes from 0.5 to 0.75 and 1.0, the non- 
robustness behavior against k becomes more and more evident. However, the slopes 
of /c/i=constant hnes remain to be smah and the increment pattern with respect to k 
is always monotonic. This suggests that the pollution error is well controlled in the 
WG solution. 




-1 -0.5 0.5 1 -1 -0.5 0.5 



Fig. 4.8. Exact solution (left) and piecewise constant WG approximation (right) for k = 100, 
and h = 1/60. 



In the rest of the paper, we shall present some numerical results for the WG 
method when applied to a challenging case of high wave numbers. In Fig. 4^ and 4.10[ 



the WG numerical solutions are plotted against the exact solution of the Helmholtz 
problem. Here we take a wave number k = 100 and mesh size h = 1/60 which is 
relatively a coarse mesh. With such a coarse mesh, the WG method can still capture 
the fast oscillation of the solution. However, the numerically predicted magnitude 
of the oscillation is slightly damped for waves away from the center when piecewise 
constant elements are employed in the WG method. Such damping can be seen in a 
trace plot along x-axis or y = 0. To see this, we consider an even worse case with 



k = 100 and h = 1/50. The result is shown in the first chart of Fig. 4.9 We 
note that the numerical solution is excellent around the center of the region, but it 
gets worse as one moves closer to the boundary. If we choose a smaller mesh size 
h = 1/120, the visual difference between the exact and WG solutions becomes very 



small, as illustrate in Fig. 4.9 If we further choose a mesh size h = 1/200, the exact 
solution and the WG approximation look very close to each other. This indicates an 
excellent convergence of the WG method when the mesh is refined. In addition to 
mesh refinement, one may also obtain a fast convergence by using high order elements 



in the WG method. Figure |4.11| illustrates a trace plot for the case of k = 100 and 
h = 1/60 when piecewise linear elements are employed in the WG method. It can be 
seen that the computational result with this relatively coarse mesh captures both the 
fast oscillation and the magnitude of the exact solution very well. 

5. Concluding Remarks. The numerical experiments indicate that the WG 
method as introduced in [54 is a very promising numerical technique for solving the 
Helmholtz equations with large wave numbers. This finite element method is robust, 
efficient, and easy to implement. On the other hand, a theoretical investigation for the 
WG method should be conducted by taking into account some useful features of the 
Helmholtz equation when special test functions are used. It would also be valuable to 
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Fig. 4.9. The trace plot along x-axis or y = form WG solution using piecewise constants. 



test the performance of the WG method when high order finite elements are employed 
to the Helmholtz equations with large wave numbers in two and three dimensional 
spaces. 
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Fig. 4.10. Exact solution (left) and piecewise linear WG approximation (right) for k = 100, 
and h = 1/60. 
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Fig. 4.11. The trace plot along x-axis or y = form WG solution using piecewise linear 
elements. 
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